Coupled Electron Ion Monte Carlo Calculations of Dense Metallic Hydrogen 



o 
o 

(N 



(N 



Carlo Pierleoni, 1 David M. Ceperley, 2 and Markus Holzmann 3 

1 INFM and Department of Physics, University of L'Aquila, Via Vetoio, 1-67010 L'Aquila, Italy 
2 University of Illinois at Urbana- Champaign, Urbana, IL 61801, USA 
:l LPTL, UMR 7600 of CNRS, Unwersite P. et M. Curie, Paris, France 

We present a new Monte Carlo method which couples Path Integral for finite temperature protons 
with Quantum Monte Carlo for ground state electrons, and we apply it to metallic hydrogen for 
pressures beyond molecular dissociation. We report data for the equation of state for temperatures 
across the melting of the proton crystal. Our data exhibit more structure and higher melting 
temperatures of the proton crystal than Car-Parrinello Molecular Dynamics results. This method 
fills the gap between high temperature electron-proton Path Integral and ground state Diffusion 
Monte Carlo methods. 
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The knowledge of the physical properties of hydrogen 
in a wide range of thermodynamic conditions is a key 
problem in planetary and high pressure physics 0,0]- In 
the search for the metallization transition three differ- 
ent insulating molecular crystal phases have been clearly 
observed so far in diamond anvil cell experiments up to 
3.2Mbar0 at room temperature and below. Metalliza- 
tion has been obtained in shock wave experiments for a 
warm dense molecular liquid 4] but properties at finite 
temperature and/or at higher pressure are largely un- 
known because experiments are increasingly difficult. 

A large body of theoretical investigations of high pres- 
sure hydrogen have appeared over the years 0. They 
helped the understanding of the experimental observa- 
tions and hold out the prospect of predicting the room 
temperature metallization pressure and the phase dia- 
gram at higher pressure. However the present under- 
standing of high pressure hydrogen is unsatisfactory be- 
cause 1) energy differences among different crystalline 
phases are small requiring a very accurate total energy 
method to determine the stable crystalline phase and lo- 
cate transition lines; 2) size effects are large in metallic 
and quasi-metallic systems and Brillouin zone sampling 
is extremely important for accurate total energy calcula- 
tions; 3) proton quantum effects are important and can 
influence the energetic ordering of crystal phases ; 4) an 
accurate theoretical prediction of metallization may re- 
quire accuracy beyond that of the LDA+GGA Density 
Functional Theory 0,0- 

Here we describe a method based on Quantum Monte 
Carlo (QMC) calculation of the electronic energy for 
quantum mechanical protons able to sample efficiently 
the protonic configurational space and spontaneously 
find the stable phase of the system within the Born- 
Oppenheimer approximation. Previous QMC studies 
of hydrogen at T = have treated electrons and pro- 
tons at the same level of description and become in- 
efficient in following the evolution of particles of very 
dissimilar mass (m p /m e — 1836). Moreover, the inter- 
esting effects of temperature are absent in this proce- 
dure. Nonetheless, they have established that pressure 
dissociation of hydrogen molecules at T=0K occurs at 



r s = [3/(4tt lie)] 1 / 3 = 1.31 (P - Mlbars) 8], where 
ri e is the electronic number density. Upon dissociation 
the molecular crystal transforms to a proton lattice of 
diamond structure and later to a lattice of cubic sym- 
metry (bec) at P > 8Mbars\ik Il0|, At finite tempera- 
ture Restricted Path Integral Monte Carlo (RPIMC) HH 
has been used to predict the equation of state (EOS) 
and to investigate the occurrence of the plasma phase 
transitionfrlj. In RPIMC, both electrons and protons 
are at finite temperature but it is efficient only for tem- 
peratures above 1/20 of the electronic Fermi tempera- 



ture (roughly 3 x 10 K at 



1). The new method 



described here, called Coupled Electronic-Ionic Monte 
Carlo (CEIMC)[lll3, is able to fill the gap between 
the RPIMC and the ground state QMC methods. We 
study metallic hydrogen in a range of densities and tem- 
peratures where molecules are absent and where protons 
undergo a solid-fluid transition. We report results for the 
EOS and give a qualitative location of the transition line. 

In the CEIMC method the proton degrees of freedom 
are advanced by a Metropolis algorithm in which the en- 
ergy difference between the actual state S and the trial 
state S' is computed by a Quantum Monte Carlo calcu- 
lation (either variational (VMC) or reptation (RQMC) 
[l^j)- The energy difference is affected by statistical 
noise which would bias the MC sampling. Unbiased sam- 
pling of the proton configurations can be achieved by the 
penalty method [lr^. a generalization of the Metropolis 
algorithm. 

We sample the electronic degrees of freedom according 
to the sum of the electronic distribution functions (e. g. 
the square of the trial wave function in VMC) for the S 
and S' states, and we compute the ener gies for the two 
states as correlated sampling averages [l3lll4| . thereby re- 
ducing the noise. Analytic trial wave functions including 
backflow and three-body correlational have been used 
in most of our calculations. These functions are particu- 
larly appropriate to our methods since: 1) they are quite 
accurate; 2) they are free of adjustable parameters so do 
not require optimization; 3) their computational cost is 
much less then solving the Kohn-Sham equations as was 
done in previous QMC calculations licj. in particular 



for a random arrangement of several tens of protons. 

To go beyond VMC, we implemented a Reptation 
Quantum Monte Carlo algorithm (RQMC)[l5j to sample 
more accurately the electronic ground state. Similar to 
Diffusion Monte Carlo (DMC), RQMC projects the trial 
wavefunction to the ground state within the Fixed-Node 
approximation. The high quality of our trial wave func- 
tions makes it possible to relax to the ground state with 
a very limited number of time slices. In RQMC the elec- 
tronic path space is sampled by a reptation algorithm in 
which, at each step, a new link is added to one end of the 
path and an existing link is deleted from the other end, 
subject to a "Metropolis" acceptance/rejection step. To 
speed up convergence we have introduced the "bounce" 
algorithm in which the gro wth direction is reversed only 
when a move is rejected|l8j. RQMC is particularly suited 
for computing energy differences. 

To reduce finite size effects in metallic systems, we av- 
erage over twisted boundary conditions (TABC) when 
computing electronic energies {i.e. we integrate over the 
Brillouin zone of the super cell) 0, ^| . All properties 
are averaged over 1000 different k-points. For the typical 
protonic displacement, we compute the energy difference 
over 100 electronic steps/k-point. After averaging over 
k-points, the noise level is small enough to simulate tem- 
peratures as low as lOOKfli). 

We represent protons by imaginary time path integrals 
without considering the statistics of the protons, (those 
effects are negligible in this temperature-density range.) 
For efficiency, it is important to minimize the number of 
protonic time slices. We have used the pair action of an 
effective proton-proton potential and treated the differ- 
ence between the true Born-Oppenheimer energy and the 
effective potential with the primitive approximation!^. 
With this action, we find that a proton imaginary time 
step t p = 0.3 x 10 _3 A" _1 is appropriate for r s > 1 so 
that few tens of time slices allow for calculations above 
100K. When coupled with TABC, we can, at each pro- 
ton move, randomly assign a subset of k-points to each 
protonic slice without introducing a detectable system- 
atic effect. This strategy allows one to simulate quan- 
tum protons at essentially the same computational cost 
as classical protons, except for the slower relaxation of 
the protonic paths. 

In order to assess the accuracy of the CEIMC method 
we first consider a system of N p — N e = 16 at 
r, = 1 and T=5000A and compare with RPIMC. 
The CEIMC-RQMC calculation, performed with r e = 
0.0125i? _1 ,/3 e = Q.5H- 1 (41 electronic time slices), pro- 
vides a total energy lower than the VMC estimate by 
A(2)mH I 'atom = 1260(630)A.However, the VMC and 
RQMC pressures agree within error bars. Comparison 
between VMC and RQMC pair correlation functions is 
also very good (see figure [TJ. The VMC and RQMC 
gep(f)'s are superimposed except at distance below 0.2ao 
(ao = 0.529 A is the Bohr radius), due to time step er- 




FIG. 1: CEIMC-RPIMC comparison for electron-proton 
and proton-proton correlation function at r, = 1,T = 
5000/^, N p = N e = 16. 



ror in RQMC. As for g pp (r), RQMC curve is slightly 
more structured than the VMC one. In RPIMC, such 
"low" temperatures (the Fermi temperature at r s = 1 
is 1.8AH = 5.8 10 5 A) can be reached only by impos- 
ing less realistic ground state nodal restriction [ill 
RPIMC data, obtained with free particle nodes and 
1000 time slices, agrees with CEIMC ones. CEIMC 
computed g ep (r) exhibits slightly more structure than 
RPIMC, and since thermal effects on the electrons should 
be largely negligible in such conditions, we attribute the 
observed difference to the more accurate nodal structure 
of CEIMC compared to RPIMC. 

Next we compare with Car-Parrinello Molecular Dy- 
namics (CPMD) simulation [2l| which uses the LDA com- 
puted forces. Figure □ shows that CEIMC- VMC g pp (r)'s 
for classical protons exhibit considerably more structure 
than does LDA. CPMD simulations considered systems 
of classical protons with a closed shell (in reciprocal 
space), and only the T point. We compare with two dif- 
ferent CEIMC calculations for classical protons, namely 
an open shell system (N p = 32) with the TABC, and a 
closed shell system (N p = 54) with the T point only. For 
the latter case, we find that the g PP (r) from VMC and 
RQMC (not shown) agree; but they exhibit more struc- 
ture than CPMD. The TABC one is in the liquid state, 
while the simulation using only the T point, initially 
prepared in a liquid state from temperature quenching, 
exhibits the onset of spontaneous crystallization. The 
larger correlation in CEIMC with respect to CPMD is 
compatible with our early estimate of the melting tem- 
perature of the fee crystal of classical proton between 
1000K and 1500K Q at variance with the the LDA esti- 
mate of 350K (for the bec crystal) [2l| . The observed dis- 
crepancy between CEIMC and CPMD is surprising since 
LDA is generally believed to be accurate at high density. 
However a previous study of hydrogen at r s = 1.31J9J 
reported that differences in energy among several crys- 
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FIG. 2: Pair correlation function at r s = 1,T = 1000A. 
Comparison between CEIMC-VMC-TABC with N p = 32, 
CEIMC-VMC-PBC N p = 54 s and CPMD-LDA N p = 162 
(simulation with N p = 54 is identical). Data from CEIMC- 
VMC-TABC at T=2000K (stars) arc also reported. 

tal structures obtained within LDA are smaller than en- 
ergy differences from Diffusion MC by roughly a factor of 
two. Also zero point energies in QMC were roughly twice 
the LDA estimates (from the harmonic approximation). 
This suggests that the Born-Oppenheimer surface from 
LDA is flatter than the more accurate one from QMC. 
Moreover there is a known issue in computing the ionic 
temperature in CPMD; the simple estimate based on the 
ionic kinetic energy provides only a lower bound for the 
true temperature |22j. Tracing the origin of the observed 
discrepancy between CEIMC and CPMD results would 
deserve an independent study. Here we just note that 
better agreement is observed between CPMD results at 
temperature T and CEIMC data at temperature 2T for 
300 <T< 3000, see for instance figure 

RQMC is roughly an order of magnitude more expen- 
sive than VMC. Therefore, it is important to establish 
the accuracy of VMC before performing a systematic 
study of the equation of state. In table [I] we compare 
VMC and RQMC at T — 5000 AT and r s = 1.2, the low- 
est density we have considered. Calculations are done at 
the r point and a projection time of (3 e = 0.68-ff -1 for 
RQMC. We have checked that for protons in bcc lattice, 
the energy has converged to its ground state value for 
(3 ~ 0.6. We find that the VMC energy is systematically 
higher by roughly 7.6(2)mH/at = 2400(60)AT,while the 
VMC pressure is systematically lower by 0.03(l)M6ars. 
The error on the energy is expected to be independent 
of the temperature and to decrease with increasing den- 
sity. Even though the amount of energy missing in VMC 
is quite large on the proton energy scale, we only ob- 
serve a minor effect on g pp (r); energy differences are 
quite accurate within VMC. On the basis of the above 
results, we performed a systematic study of the Equa- 
tion of State (EOS) using VMC. In table HJ we report 



TABLE I: Comparison between VMC and RQMC energies 
and pressures for r s = 1.2, T = 5000A, N p = N e = 54, PBC. 
r = is the estimate with RQMC time step errors removed 
by extrapolation, a 2 is the variance of the local energy in 
VMC. 



Te 


E fot (h/at) a 1 E kin E pot P (Mbars) 


vmc 
0.01 
0.00 


-0.4694(2) 0.0472(4) 0.8812(4) -1.3508(4) 5.55(1) 
-0.4768(4) 0.8850(6) -1.3618(6) 5.50(1) 
-0.47696 0.89112 -1.36808 5.581 



total, kinetic and potential energies, pressure, the Lin- 
demann ratio for bcc crystal, and the proton kinetic en- 
ergy. The latter quantity can be compared to 3ATbT/2 
(last column). The zero point proton motion affects not 
only the proton kinetic energy but also increases the elec- 
tronic kinetic energy and, to a smaller extent, the con- 
figurational energy. At r s = 1 and T — BOOK we find 
a total energy increase of \A.%(2)mH / at — 4670(60) AT of 
which 2020(30)Af comes from the proton kinetic energy, 
2200(20)AT the electronic kinetic energy, and 450(10)^ 
the configurational energy. Residual finite size effects 
have been estimated from static lattice calculation at 
r s = 1 to be of the order of lOmH/at on the energy, 
and 0.21 Mbars on the pressure. The transition line, esti- 
mated by the dynamical observation of melting, is located 
between 1000K and 2000K at r s = 0.8, between 500K 
and 1000K at r s = 1.0 and close to 1000K at r s = 1.2. 
Indeed at the latter density and at T=1000K the system 
is able to sustain both liquid and crystal states for the 
entire length of our simulations (80000 protonic steps). 

In conclusion we have developed a new and efficient 
Quantum Monte Carlo Method to study low temperature 
quantum protons and ground state electrons which is a 
major improvement over previous QMC and DFT-LDA 
based methods. It allows for simulations of many-body 
hydrogen using QMC for the electronic energies. We have 
developed efficient procedures to include protonic path 
integrals and k-point sampling. We have applied it to 
metallic hydrogen beyond molecular dissociation and in- 
vestigated the solid- fluid transition of the protons. The 
present methodology can be extended in several ways. 
Constant pressure algorithm would be useful to study 
structural phase transitions. However for metallic sys- 
tems, we have found that level crossings, arising from 
changes in the shape of the simulation box, considerably 
increase the noise level and makes our correlated sam- 
pling procedure inefficient. The method can be easily 
extended to the insulating molecular phase by replac- 
ing the metallic trial functions with localized molecular 
orbitals[l3l flij]. A study of the melting line of molecu- 
lar hydrogen is in progress. Consideration of the metal- 
insulator transition requires a trial function that goes 
smoothly from metallic to localized orbitals. We are in- 
vestigating an accurate and efficient form for this. Ex- 
tension of the present method to more complex elements 
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TABLE II: Energy and pressure for a system of N p = 54 quantum protons with VMC-TABC. Units of energy are 
hartrees/proton. M v is the number of protonic time slices (M p = 1 means classical protons). 7 is the rms deviation di- 
vided by the nearest neighbor distance for a bcc lattice. 





T(KK) 


M p 


E 




Epot 


P(Mbars) 




K p x 10 2 


Kp x 10 2 


0.8 


0.5 


16 


-0.0594(2) 


1.8419(1) 


-1.9033(1) 


81.07(3) 


0.169(1) 


1.57(3) 


0.2375 




1.0 


16 


-0.0586(4) 


1.8428(4) 


-1.9034(1) 


81.16(3) 


0.183(1) 


1.53(4) 


0.475 




2.0 


8 


-0.0522(4) 


1.8338(4) 


-1.9018(1) 


81.69(3) 





1.78(3) 


0.950 




3.0 


4 


-0.0442(4) 


1.8538(6) 


-1.9000(2) 


82.33(6) 





2.14(7) 


1.425 




4.0 


4 


-0.0382(8) 


1.8590(8) 


-1.8991(1) 


82.83(6) 





2.57(7) 


1.900 




6.0 


2 


-0.0268(8) 


1.8688(8) 


-1.8974(2) 


83.80(6) 





3.29(4) 


2.850 




10.0 


1 


0.016(1) 


1.8886(8) 


-1.8934(4) 


85.78(9) 





4.750 


4.750 


1.0 


0.5 


8 


-0.3512(2) 


1.2142(2) 


-1.5655(1) 


20.101(3) 


0.177(1) 


0.97(2) 


0.2375 




1.0 


4 


-0.3480(2) 


1.2176(2) 


-1.5657(1) 


19.68(1) 




1.07(2) 


0.475 




2.0 


4 


-0.3430(2) 


1.2260(4) 


-1.5653(1) 


20.65(1) 





1.44(2) 


0.950 




3.0 


2 


-0.3356(4) 


1.2298(4) 


-1.5655(1) 


20.83(1) 




1.72(3) 


1.425 




5.0 


1 


-0.3262(6) 


1.2390(6) 


-1.5652(1) 


21.26(2) 




2.375 


2.375 




10.0 


1 


-0.2888(6) 


1.2740(4) 


-1.5630(2) 


22.95(3) 




4.750 


4.750 


1.2 


0.3 


10 


-0.46610(4) 0.8776(1) 


-1.3437(1) 


5.554(1) 


0.134(1) 


0.59(1) 


0.1425 




0.5 


8 


-0.4661(1) 


0.8792(1) 


-1.3439(1) 


5.594(3) 


0.177(2) 


0.67(1) 


0.2375 




1.0 


4 


-0.4632(1) 


0.8811(2) 


-1.3443(2) 


5.641(3) 


0.196(3) 


0.77(1) 


0.475 




1.0 


4 


-0.4610(2) 


0.8858(2) 


-1.3468(1) 


5.735(6) 


liquid 


0.77(1) 


0.475 




2.0 


4 


-0.4552(2) 


0.8918(2) 


-1.3469(1) 


5.893(6) 




1.19(3) 


0.950 




3.0 


2 


-0.4492(4) 


0.8996(3) 


-1.3488(1) 


6.08(2) 




1.53(3) 


1.425 




5.0 


1 


-0.4386(6) 


0.9106(4) 


-1.3492(2) 


6.37(2) 




2.375 


2.375 




10.0 


1 


-0.4036(6) 


0.9478(4) 


-1.3514(1) 


7.34(2) 




4.750 


4.750 



is straightforward, provided we have efficient trial func- 
tions. 
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